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Abstract 

In high energy physics, a widely used method to treat systematic uncertainties in 
confidence interval calculations is based on combining a frequentist construction of 
confidence belts with a Bayesian treatment of systematic uncertainties. In this note 
we present a study of the coverage of this method for the standard Likelihood Ratio 
(aka Feldman & Cousins) construction for a Poisson process with known background 
and Gaussian or log-Normal distributed uncertainties in the background or signal 
efficiency. For uncertainties in the signal efficiency of upto 40 % we find over-coverage 
on the level of 2 to 4 % depending on the size of uncertainties and the region in signal 
space. Uncertainties in the background generally have smaller effect on the coverage. 
A considerable smoothing of the coverage curves is observed. A software package 
is presented which allows fast calculation of the confidence intervals for a variety 
of assumptions on shape and size of systematic uncertainties for different nuisance 
parameters. The calculation speed allows experimenters to test the coverage for 
their specific conditions. 
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1 Introduction 

The calculation of confidence intervals in case of presence of systematic uncer- 
tainties is an open problem. Systematic uncertainties are those uncertainties 
present in parameters which will affect the calculated confidence interval, but 
which are not of prime interest, so called nuisance parameters. Examples for 
nuisance parameters are the efficiency or the predicted background rate. 
In 1992 Cousins & Highland [1] proposed a method which is based on a 
Bayesian treatment of the nuisance parameters. They proposed to perform 
a frequentist construction (following Neyman [2] ) and to replace the probabil- 
ity density function (PDF) describing the statistical process of prime interest 
with a PDF which is obtained by a convolution the original PDF with the one 
describing the uncertainties in the nuisance parameters: 



where e' is the true value of the nuisance parameter, e denotes its estimate 
and s and n symbolize the signal hypothesis and the experimental outcome 
respectively. 

Highland & Cousins only treated the case of Gaussian uncertainties in the 

signal efficency. The method has since been generahzed to operate with the 
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modern unified ordering scheme proposed by Feldman & Cousins [3] and tak- 
ing into account several nuisance parameters (efficiencies and background) and 
correlations [5]. This generalized method has already been used in a number 
of particle and astroparticle physics experiments (e.g. [6] [7] [8], [9], [10]). 
The most crucial property of methods for confidence interval construction is 
the coverage, which states: 

A method is said to have coverage (1-a) if, in infinitely many repeated experi- 
ments the resulting confidence interval includes the true value with probability 
(1-a) irrespective of what the true value is. 

100(1-0;)% is hereby commonly taken to be 68% 90%, 95% etc. Recently, 
the coverage properties of a fully frcquentist method, the Profile Likelihood 
method, have been studied [11]. The Profile Likelihood method was found to 
have surprisingly good coverage properties with a small (mostly negligible) 

amount of under cover age. 

In this note we undertake a systematic study of the coverage of the Bayesian 
method in order to give more substance to further recommendations on what 
method to use to calculate confidence intervals in presence of systematic uncer- 
tainties. Previous studies [12] of this method dealt only with certain limiting 
cases and were constrained by computational requirements. 
The note is organized as follows: in the next section we will review the frc- 
quentist construction of confidence intervals, in particular the likelihood ratio 
ordering scheme. We will describe the Bayesian method to incorporate sys- 
tematics in section 3. The C-|— I- library used to perform the calculations is 
described in section 4. Coverage tests and connected subtleties are discussed 
in section 5. The final section is devoted to discussion and conclusion. 
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2 Frequentist confidence interval construction 



Let us consider a Poisson probability density function (PDF), p{n)s+b, for a 
fixed but unknown signal, s, in the presence of a known background with mean 
b. For every value of s we can find two values rii and n2 such that 



where 1 — a denotes the confidence level (usually quoted as a 100(1-q;)% 
confidence interval). Since we assume a Poisson distribution, the equality will 
generally not be fulfilled exactly A set of intervals [ni (s + b,a), n2 {s + b,a)] is 
called a confidence belt. Graphically, upon a measurement, Uo, the confidence 
interval [si,S2] is determined by the intersection of the vertical line drawn 
from the measured value Uo and the boundary of the confidence belt. The 
probability that the confidence interval will contain the true value s is 1 — a, 
since this is true for all s per construction. The choice of the rii and 712 is, 
however, not unique to define the confidence belt. An additional criterion 
has to be apphed. The currently recommended ordering scheme [3] [4] of 
the elements in the sum in equation (2) is based on likelihood ratios. This 
approach automatically provides central confidence intervals when motivated 
and upper limits when necessary, therefore it is often denoted as the "unified 
approach". The following algorithm is applied in solving equation (2): 
For each n the Sbest is found which maximizes the likelihood jC{n)s+b- In case 
of a simple Poisson distribution with known background, Sbest is given by 
■max{0, n — b). Then for a fixed s the ratio 



p{n')s+b = 1-0; 



(2) 



n'=n\ 



R{s,n)c 




(3) 
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is computed for each n, and all n's are consequently ranked according to the 
value of this ratio. Values of n are included in the confidence belt starting 
with the n with the highest rank (largest Rc) and then decreasing rank until 
J2n=ni = 1 — a. After the confidence belt has been constructed in this 

way, the confidence interval [si, S2] is found as described above. 



3 Incorporation of systematic uncertainties: the Bayesian way 



The ordering schemes are unaffected by the way of treating systematic uncer- 
tainties considered here. As mentioned earlier, the PDF describing the statisti- 
cal process will however be modified. Two concrete examples are the following: 
In the case that the only uncertainty present is a theoretical (assumed Gauss- 
shaped) uncertainty of the background process the PDF is modified to: 



00 



q{n)s+b = -j=- J p{n)s+b' e dh' (4) 



Here h is the estimated background level, and db is the uncertainty in the 
background estimation. If, in addition to the theoretical uncertainty for back- 
ground, there is the need to include the uncertainty in the signal detection 
efficiency the expression for q{n)s-[.h might be extended to: 



(5) 



/o°°/o°°P(^)6'+.'. e e dh'de' 



where 0"^ is the uncertainty in the detection efficiency expressed in relative 

terms with respect to the nominal efficiency ^ . It is important to realize that 
^ of course the eflaciency could be defined in absolute terms as well. 
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the integration variables, here e' and 6', are the possible "true" (but unknown) 
values of nuisance parameter. This indicates that this method is based on 
Bayesian statistics. 

Some examples for confidence intervals computed by this method are shown 
in table 1. 

4 pole++ 

For the coverage studies presented in this paper a reasonably fast and efficient 
code is required. Hence, a user-friendly and flexible C-|— I- library of classes 
was developed based on the FORTRAN routine presented in [13]. The hbrary 
is independent of external hbraries and consists of two main classes, Pole and 
Coverage. The first class takes as input the number of observed events, the ef- 
ficiency and background with uncertainties and calculates the limits using the 
method described in this paper. The integrals are solved analytically. Cov- 
erage generates user-defined pseudoexperiments and calculates the coverage 
using Pole. Presently the library supports Gauss and log-Normal distributed 
PDF for description of the nuisance parameters. Flat and used-defined PDFs 
are about to be implemented as well as correlations for the Gauss case. The 
class is dynamically optimized depending on if one wishes to calculate sin- 
gle (or few) confidence intervals or if one wants to perform a coverage study. 
Without these optimisations the calculation of a single interval takes about 1 
second (wall clock time) on a 1 GHz Pentium HI processor. The duration of 
a full coverage study (typically requiring the calculation of 0(10^ ) confidence 
intervals) ranges between a couple of minutes for small uncertainties and small 
signal hypotheses to order of 10 hours for large uncertainties and high signal 
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hypotheses. The perfomance of a coverage study thus seems feasable for the 
particular set of systematic uncertainties that may appear in real experiments. 
The pole++ library can be obtained from http: / /cern.ch/tegen/ statistics. html 



5 Coverage Studies 

The coverage of the method is studied using MC simulations of pseudo- 
experiments with given true value of the prime parameter (signal) or nuisance 
parameter (efficiency, background). The estimated values of the nuisance pa- 
rameters were assumed to be Gaussian or log-Normal distributed around the 
given true value. The outcome of one experiment thus consisted of a number 
of observed events (following a Poisson distribution with known background 
and depending on the true efficiency and true background) and the estimate 
of the nuisance parameter (s). 

Figure 1 shows the coverage^ as a function of signal hypothesis for two dif- 
ferent sizes of uncertainties (5 % and 40 %). The uncertainty considered is in 
the signal efficiency and assumed to be Gauss-shaped. It can be seen that the 
Bayesian method causes over-coverage which is larger for larger systematic 
uncertainty. This is also reflected in figure 3 where the mean coverage (mean 
taken over all tested signal hypotheses) is shown as a function of assumed 
uncertainty. For Gaussian uncertainties in the signal efficiency we find an in- 
crease in mean coverage by ~ 1% in the uncertainty range between 10 % and 
40 %. 

The see-saw structure which is generally seen in the no-uncertainty case (due 

^ we will denote the coverage calculated by the MC simulations as (1 — a)e// to 
distinguish it from the nominal coverage (1 — a) 
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to the discrete experimental outcome) is considerably smoothened when sys- 
tematic uncertainties are introduced. This is exemplified in figure 2, which 
compares the coverage curve with 5 % uncertainties on the signal efficiency 
with the zero uncertainties case. To quantify this further we present the rms 
of the coverage together with the mean in figure 3. The smoothing is due to 
the fact that we add a continuous variable to the problem, meaning there are 
more degrees of freedom to fulfil the sum condition of equation 2. The effect 
of the smoothing is that whereas for some signal hypothesis the coverage is in- 
creased for others it is decreased with respect to the zero uncertainty case. The 
mean coverage is therefore only rather weakly dependent on the systematic 
uncertainties. For higher signal hypotheses, where the effect of the smoothing 
is less pronounced, the effect of the uncertainties is therefore stronger. Taking 
the mean over signal hypotheses s^j-^g >6 the mean coverage increases from 92 
% (at zero uncertainties) to 94 % (at 40 % uncertainties). A side effect is that 
the introduction of rather small uncertainties seems to improve the coverage 
with respect to the zero uncertainty case for parts of the tested hypothesis 
space (see e.g. figure 2). 

Figure 4 shows the coverage for two different sizes of Gaussian uncertainty 
on the background estimate. Except for the smoothing effect the coverage is 
seemingly independent of the size of the uncertainties. An uncertainty in the 
background will yield similar results to an uncertainty in the signal efficiency 
only in the regime were the signal hypotheses are of similar size as the back- 
ground expectation. For larger signal hypotheses the coverage curve for the 
uncertainty in the background will be approaching the zero uncertainty case, 
which is refiected in a slight slope of the coverage curve for high uncertainty 
in figure 4. The result is thus that a mean coverage depending on uncertain- 
ties in the background will be less affected than the corresponding curve with 
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uncertainties in the signal efficiency. We sliow tlie mean coverage as well as 
the rms in figure 5. 

In figure 6 the effect of having to consider uncertainties both in signal and 
background detection efficiency is visible. At large signal hypotheses the cov- 
erage plot is dominated by the uncertainty in signal efficiency, at low signal 
hypotheses coverage benefits from the smoothing effect, since we added yet 
another degree of freedom. 

It should be noted that for large uncertainties the Gaussian model is not ap- 
propriate. The reason is that a Gaussian model will then significantly extent to 
the unphysical region of the space of experimental outcomes. In the Bayesian 
treatment this case is dealt with by truncating the Gauss distribution at zero 
^ , and this is consequently the way it is dealt with during the coverage test. 
Considering the coverage test there is an additional subtlety: in order for the 
measured efficiency to be a maximum likelihood estimate of the true efficiency 
we would have to renormalize the Gauss distribution to the truncated Gauss 
distribution in order to obtain the correct PDF. However, instead of doing this 
a posteriori fix, it is more reasonable to use a log-Normal distribution to model 
the uncertainty in the efficiency. The coverage for the log-Normal distribution 
is shown in figure 7. As can be seen for the highest uncertainties considered 
in this note, the Gauss distribution is still a very good approximation to the 
log Normal model. 



^ it is worth noting that Cousins &; Highland develope their method for unbiased 
estimators, i.e. the problem would not be treated in this way. 
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6 Discussion and Conclusions 



In this note we presented coverage studies for the Bayesian treatment of sys- 
tematic uncertainties. One overall conclusion is that the Bayesian treatment 
leads some over-coverage. However, introducing a continuous nuisance param- 
eter into the discrete Poisson problem results in a smoothing of the coverage 
curves. The mean coverage is therefore only weakly affected by the Bayesian 
treatment of nuisance parameters and under certain circumstances even im- 
proves with respect to the zero uncertainty case. 

In a frequentist approach the meaured estimate of a nuisance parameter is 
considered to be distributed around a given true value, which is consequently 
the way coverage was calculated. The Bayesian method on the other hand 
views the true value as distributed around the measured value. The underly- 
ing assumption for going from one approach to the other is (at least in case 
of a Gauss-distribution) a flat prior probability of hypotheses. The present 
study indicates that this assumption does not lead to a violation (except for 
over-coverage) of the coverage requirement. 

The routines used for the presented calculations are reasonably fast and pub- 
licly available. A coverage study is therefore feasable for each problem at hand. 
The confidence level required for the confidence intervals could then in prin- 
ciple be adjusted to recover correct coverage. 
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8 Figures and tables 
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Table 1 

Examples of likelihood ratio 90% confidence intervals with Bayesian treatment of 
systematic uncertainties. Uncertainties are assumed to be Gauss distributed in the 
signal efficiency. 
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Fig. 1. Calculated coverage as function of signal hypothesis. Two case are shown: 5 
% and 40 % Gaussian uncertainties in the signal efficiency. The nominal coverage 
was 90%. 
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Fig. 2. Calculated coverage with 5 % systematic uncertainties is compared with the 
zero uncertainties case. The nominal coverage was 90 %. 
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Fig. 3. Calculated mean coverage and rms as function of Gaussian shape uncertain- 
ties in the signal efficiency. The mean was calculated over the full range of signal 
hypotheses and for signal hypotheses larger than 6, respectively. The rms was cal- 
culated over the full range of signal hypotheses. The nominal coverage was 90 %. 
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Fig. 4. Calculated coverage as function of signal hypothesis. Two case are shown: 
5 % and 40 % Gaussian uncertainties in the background prediction. The nominal 
coverage was 90 %. 
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Fig. 5. Calculated mean coverage and rms as function of Gaussian shape uncertain- 
ties in the background prediction. The mean was calculated over the full range of 
signal hypotheses and for signal hypotheses larger than 6, respectively. The rms was 
calculated over the full range of signal hypotheses. The nominal coverage was 90 % 
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Fig. 6. Coverage for the case that uncertainties are present in both the signal de- 
tection efficiency and in the background prediction. The nominal coverage was 90 
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Fig. 7. Coverage as a function of signal hypothesis for a Gaussian PDF and a log 
Normal PDF 
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